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Abstract 

We give a constructive method for realizing an arbitrary directed graph (with no 
one-cycles) as a heteroclinic or an excitable dynamic network in the phase space of a 
system of coupled cells of two types. In each case, the system is expressed as a system 
of first order differential equations. One of the cell types (the p-cells) interacts by 
mutual inhibition and classifies which vertex (state) we are currently close to, while 
the other cell type (the y-cells) excites the p-cells selectively and becomes active only 
when there is a transition between vertices. We exhibit open sets of parameter values 
such that these dynamical networks exist and demonstrate via numerical simulation 
that they can be attractors for suitably chosen parameters. 


1 Introduction 

Researchers in neuroscience often want to understand how the structure of the central ner¬ 
vous system of animals relates to the function of the system both in healthy and diseased 
individuals, and much effort has been spent trying to model coupled neurons as networks 
of nonlinearly interacting cells (see e.g. [i5i[in]) and emergent dynamical properties of the 
network are clearly important for an understanding of neural function from basic signal 
processing to high-level cognition; for example m suggest that the “Dynome” of possible 
states of the network is just as critical as the physical “Connectome” of connections between 
neurons and neural assemblies. 
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In addition to the physical network of conpled cells, the possible transient dynamical 
states of the system may usefnlly be thonght to form a network in phase space of some type. 
A heteroclinic network is a hnite set of dynamical states connected by trajectories (dehned 
more precisely below). This idea has been developed in a number of models inspired 
by neural systems HI El 1271 [19]; see also winnerless competition dynamics |2Tl |5], stable 
heteroclinic channels jH], and networks of unstable attractors [23] . If there are no “direct 
connections” between states, there may be excitable connections where a certain amplitude 
of a on-off perturbation at a state is needed for a transition from one state to another, giving 
rise to an excitable network (also dehned below) whose structure will typically depend on 
the amplitude of perturbation usedj^ 

In a previous paper [B], we proposed two constructions to show how arbitrary graphs 
may be embedded or realised as a heteroclinic network in the sense that there is a one-to-one 
mapping between vertices and edges of the graph and dynamical states. That paper presents 
a “simplex network” and a “cylinder network” of coupled cells that allows one to realise any 
hnite directed graph that is one-cycle free into phase space as a heteroclinic network. The 
vertices of the graph correspond to equilibria that are saddles in phase space and the edges of 
the graph correspond to connecting (heteroclinic) orbits in phase space. The constructions 
in [6] require several diherent cell types; Field na has recently shown that it is possible to 
do this even if one restricts to just one cell type. 

This paper aims to present an explicit construction to realize arbitrary graphs in phase 
space as either heteroclinic or excitable networks. Not only this, we identify a bifurcation 
from a heteroclinic network to an excitable network of the same topology on changing a single 
parameter in the governing equations. In the former case the vertices are saddle equilibria 
and the connections are heteroclinic. In the latter case the vertices correspond to stable 
equilibria in the network that are sensitive to perturbations in directions corresponding to 
the edges in the graph. The network construction uses two cell types, where there is strong 
inhibition between cells of one type and strong excitation of this cell type by the second cell 
type, reminiscent of neuronal systems. 

We structure the paper as follows: after dehning heteroclinic and excitable networks in 
phase space we give a simple motivating example (given in equation ([^) of a system that 
realises a cyclic graph between three nodes as either a heteroclinic or an excitable network, 
depending on values of the parameters in the equations. In Section]^ we introduce an explicit 
description of a coupled cell model (given in equation ([^) that is parameterised by a number 
of constants. The hrst main result. Theorem 2^, shows that an arbitrary hnite directed 
network can be robustly realised as a heteroclinic network in the phase space of this coupled 
cell system. Similarly, the second main result. Theorem |2.4[ shows that an arbitrary hnite 


directed network can be robustly realised as an excitable network for amplitude 6 in the 
phase space of this coupled cell system. The minimum amplitude is related to the distance 
of a parameter from a bifurcation point where the saddles in the heteroclinic network are 


^Excitable networks in this sense have been considered before (for example see (Sj), but should not be 
confused with networks of coupled excitable units (for example see |22jl. that may or may not have excitable 
networks in phase space, depending on the coupling. 
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stabilized. In particular, the 6 may be made as small as desired by choosing parameters 
appropriately. 

We give some numerical examples in Section that realise the Kirk-Silber network [IH] 
of competition between two cycles using (|^ in the two different ways outlined in Section 
In the presence of noise, we note that typical trajectories explore either cycle in a random 
manner. In Section we discuss some implications of the study, including generalizations 
that give networks in phase space where there may be a mixture of heteroclinic and excitable 
connections, and where the excitable connections may have various thresholds. 

1.1 Heteroclinic and excitable networks in phase space 

Consider an ODE with phase space x G M'^, dehned by 

^ = x = /(x), (1) 

and suppose that the flow generated by the solution of this ODE starting at xq is x{t) = 
0t(xo). Let Bs{x) denote the closed ball centered on x with radius 5 > 0. For ^ an equilibrium 
of Q we dehne the stable and unstable sets lE^(^) = {y : |0i(|/) — 0 as t —)■ oo} and 

iy“(^) = {y : \(l>t{y) — ^1 —)■ 0 as t —)• —C)o}; these are manifolds if ^ is hyperbolic. Typically 
we will consider the case that all equilibria are hyperbolic. 

Our dehnition of a heteroclinic network is substantially weaker than that given in most 
of the literature (see below for further details); we say a set X C is a (weak) heteroclinic 
network (in phase spacej|^ if there is a set of equilibria such that 

n 

X = Xhet({e.}) := U n 

*J=1 

and we say there is a heteroclinic connection from to whenever 

We assume there are no homoclinic connections, i.e. that nlE®(.^j) = {■Ci}- Note that 

this dehnition of heteroclinic network is weaker than that used in most of the literature, 
e.g. dH IB], in the following ways: (a) we do not require any chain recurrence or even 
connectedness of the network; for example, we do not exclude the possibility that the system 
is of gradient type; (b) we do not require that the entire unstable set is contained in the 
network; (c) we do not require that the equilibria are hyperbolic, although in typical cases the 
equilibria of heteroclinic networks are saddles (if there are incoming and outgoing heteroclinic 
connections at that equilibria), and the equilibria of proper excitable networks are sinks. 

^We refer to a “(weak) heteroclinic network (in phase space)” simply as a “heteroclinic network” for the 
remainder of the paper. 
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We say a set X C is an excitable network (in phase space) for amplitude 5 > 0 if there 
is a set of equilibria such that 

n 

X = Xe.c({6}, := U ^ ^ ^ > 0 } n 

*J=1 

In other words, an excitable network is the union of a number of equilibria and the set 
of trajectories within the stable manifolds of these equilibria that come within S of other 
equilibria. We say there is an excitable connection for amplitude 5 > 0 from fi to fj whenever 

We say an excitable connection from fi to fj has threshold Sth{^i,^j) where 

:= inf{<5 > 0 : n ^ 0}. (2) 


An excitable network for amplitude 6 is proper if all of its excitable connections have hnite 
threshold, i.e. if there is a 5', with 5 > 5' > 0 such that there are no excitable connections 
for amplitude 5' from any fi to any fj within the network. In this terminology, a hetero¬ 
clinic connection from fi to corresponds to there being an excitable connection with zero 
threshold. 

An excitable network X is forwards (but not necessarily backwards) invariant {(j)t{X) C 
X for all f > 0) while a heteroclinic network is both forwards and backwards invariant 
{(j)t{X) C X for all f G M). Observe also that any hnite set of equilibria are connected in an 
all-to-all manner by taking excitable connections with a large enough amplitude. 

Now consider a (hnite) graph T = (V,£^) with vertices V = {t)i,..., and 

directed edges £ = {ci,..., e„^}. We dehne a{k) and oj{k) so that is the edge from t)a{k) 
to We say T is one-cycle free if a{k) ^ u{k) for all k and we will assume henceforth 

that r is one-cycle free. We say a heteroclinic network X realises the graph T if each vertex 
Dj of r corresponds to an equilibrium in X, and there is an edge of T from Dj to Dj if and 
only if there is a connection from fi to fj in X. We say an excitable network X for amplitude 
5 realises the graph F if each vertex Dj in F corresponds to an equilibrium in X and there 
is an edge in F from to dj if and only if there is a connection in X for amplitude 6 from 
to fj. In Section]^ theorems 2.3 and 2.4, we present an explicit system whereby any graph 
F can be realized as a heteroclinic network or as an excitable network for some small 5 > 0. 

Note that for a given ODE and set of equilibria there can be a mixture of heteroclinic 
and excitable connections with diherent thresholds. Figure [^illustrates this and shows how 
choosing different amplitudes may give rise to excitable networks with differing topology. 

There is a subtle difference between existence of a heteroclinic connection between two 
equilibria and existence of an excitable connection with threshold zero. More precisely one 
can show the following difference: 


Lemma 1.1 Consider a given ODE and two equilibria ^ 2 - 
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If we examine the three equilibria ^2 and ^3 there are heteroclinic connections from ^2 to ^3 
and from ^3 to ^ 1 . There is an excitable connection with threshold Sth > 0 that corresponds 
to the radius of the black circle around ^ 1 ; for any 6 > 6th (see for example the red circle) 
there is a connection shown in red from ^1 to ^2 with amplitude 5. (b) and (c) show the 
excitable networks between the for amplitudes 5 < Sth and S > Sth respectively; note 
that the heteroclinic network between these equilibria is (b). The threshold corresponds to 
the distance of the stable manifold (shown by the dotted line in (a)), of the saddle equilibrium 
C, from ^ 1 . 


• There is a heteroclinic connection from to f 2 if and only if there is a trajectory x(t) 
such that 


= lim x(t), and ^2 = lini x(t). 

t^ — OO t^+OO 


• There is an excitable connection from ,^1 to f 2 with threshold zero if and only if there 
are trajectories xsit), S > 0 such that 


lim lim inf \xs{t) — ^i| = 0 , and ^,2 = lim xsit) for all 5 > 0 . 

5^0 t^ — OO t^+cxo 


Proof: This follow from considering the dehnitions of heteroclinic and excitable connection. 
If there is an excitable connection then for arbitrarily small amplitudes S there for each S 
we have a trajectory xsit) that approaches closely to .^1 in the limit t —>■ — 00 ; however one 
may need to choose a different trajectory on reducing S. QED 

Figure [^illustrates that an excitable connection with threshold zero may be a connection 
of “depth two” or greater [2] even if one can take the same trajectory independent of S in 
Lemma 11.11 


1.2 Example: a cycle of order three 

We now give a motivating example of how a directed graph can be used to design a coupled 
cell system. Consider the three-node, three-edge cyclic graph in Figure [^a). Using two 
types of dynamical cells we construct a system consisting of six cells (given by equation (j^), 
as shown in Figure [^b), where the full coupling between the cells is shown in Figure [Sp^). 
The p-cells classify the location (one p cell is active at each vertex in the graph), while the 


5 






Figure 2: Schematic diagram illustrating the difference between a heteroclinic connection 
and an excitable connection with zero threshold. Four equilibria are such that there is an 
excitable (but not a heteroclinic) connection shown in red from any of the (^j, z = 1, 2,3 to .^4. 
Note that the alpha-limit set of the red trajectory contains the heteroclinic cycle between 
the i = 1,2,3; it is a “depth two” connection [2]. 


?/-cells only become active during transition between vertices. 

Pi = + Diplp"^ - /)) + E{-ylpiP2 + vIpI) + rjpWi 

P2 = P2{F{1 - p^) + D{pIp^ - /)) + E{-ylp2P3 + ylpl) + PpW2 

Ps = P3{F{1 - p^) + D{pIp^ - p^)) + E{-ylp^pi + ylpl) + ppW^ 

yi = 9{yi:A - Bpj + C{y^ - yl)) + yyW^ 
m = p(l/2. A- Bpl + C{y^ - yl)) + PyW^ 
m = g{y3, A-Bpl + C{y^ - yl)) + PyWe 


The Wj are white noise processes, Pp and py are noise amplitudes and the function g is given 
in equation ([^. We choose a standard set of parameters (these lie within an open region 
(12) of suitable parameters described in section]^, and consider the effect of low amplitude 
noise: 


A = 0.5, B = 1.8, C = 2, D = 10, E = 4, E = 2, pp = py = 10 


-3 


(4) 


The connections between the cells in ® are mostly inhibitory (negative feedback), except 
for the connections shown in Figure l^c) as solid which represent excitatory connections 
between selected cells in the sense that they provide positive feedback. Theorem 2.4 from 


the next section can be used to deduce that there is a heteroclinic cycle as shown schemat¬ 
ically in Figure [^d). Finally, for the same parameters as in (|^ except choosing B = 1.49, 
Theorem |2.3| shows that there is an excitable network as shown schematically in Figure [^e). 

Figure illustrates the attracting behaviour of this system: in the absence of noise 
(a) and (c), the behaviour of the heteroclinic and excitable networks are quite different. 
In the presence of noise (b) and (d), they are qualitatively similar due to the trajectories 
being driven around the network by the noise. In hgure we show detailed time-series of 
the system, illustrating the transitions corresponding to edges between the vertices of the 
directed graph. Note that yi is switched on during the transition from (pi = 1) to ^2 
{p 2 = 1). Throughout this paper, we use a Heun integrator with timestep h = 0.01 for 
simulations of the noise-driven systems. 
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Figure 3: (a) The cyclic graph: three edges e* connect three vertices dp (b) schematic coupling 
architecture of the six-cell network realizing (a). The p-cells classify the location when at 
the vertices of (a) while the |/-cells only become active during transition between vertices, 
(c) shows all the connections between the cells; the dashed arrows indicate inhibitory while 
the solid arrows indicate excitatory connections. (d,e) schematically show the connections in 
phase space for this network, where in (d) are saddles connected by heteroclinic connections 
and (e) are stable nodes that are connected by excitable connections for amplitude 6, with 
the separatrices being the stable manifolds of the saddles Q that are close to the 


2 The coupled system with two cell types 


Consider a directed graph T = (V,T), with vertices and edges dehned as in section 1.2 


We will realise this both as heteroclinic and as excitable networks in the phase space of a 
set of coupled cells of two types: the p-cells are associated with the vertices V while the 
y-ceWs are associated with the edges £. The system we consider has phase space and 

coordinates (p, y) = (pi,... ,Pn„, Pi, • • •, Pne) ^ governed by: 


= pj(F(l -p2) + D{p]p^ -p"^)) + E{-Zf\p,y) + Zf{p,y)) 
j^Vk = 9 {yk, A - - yl)) 


for j = 1, • ■ ■ , and /c = 1, ■ ■ ■ , rig, where p^ = Pp / = EJli P% P^ 
A, B, C, D, E are constants. The function g is dehned by 


E Ue 

.i=iy. 


^ and 


giVk, A) = -yt ((St - 1)^ + A) 


( 6 ) 


while the inputs to the pj cells from the p cells are: 

zf\p,y)= -ylPc.ik)Pj 

{k : aik)=j} 

Z^\p,y)= yl'Pl{.k'y 

{k' : 

Equations ([^ have equilibria at = (0,..., 1,..., 0) G for j = 1,..., where the 

“1” is in the jth position. That is, the equilibria are at points corresponding to unit vectors 
where one of the pj is non-zero. 
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Figure 4: Trajectories for the coupled cell system ([^, for four different parameter sets. In 
each panel, pi is shown by a blue solid line, p 2 by a red dashed line, and ps by a black dotted 
line. The yj components are not shown: see hgure (a) shows the trajectory approaching 
a heteroclinic network; B = 2.5 and no noise; Vp = % = 0. Note the cycling between 
three states while slowing down typical of a heteroclinic cycle attractor, (b) is as (a) except 
with non-zero noise; Vp = Vy = 10“^. Observe that the slowing down is replaced by an 
approximate periodicity induced by the noise. In (c), parameters are chosen so there exists 
an excitable network with no noise; B = 1.49, Pp = Py = 0. The trajectory approaches a 
stable equilibrium that depends on initial conditions, (d) is as in (c) except for non-zero 
noise; Pp = 10“^, py = 5 x 10“^. Here, the noise pushes the trajectory over the thresholds 
at each equilibria and cycling behaviour is seen. (Parameters are as in 0 , except where 
stated.) 
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Figure 5: Time series of transitions between and ^2- In each panel, pi is shown as a blue 
solid line, p 2 as red dashed and yi as black dotted, (a) has B = 1 . 8 , rjp = Py = 10 “^, and 
corresponds to a heteroclinic connection, (b) has B = 1 . 49 , Pp = 10 “^, py = 3 x 10 “^ and 
corresponds to an excitable connection. (Parameters are as in 0 , except where stated) 

Note that y = g{y, A) has a hysteresis loop that can be switched by changing A through 
the interval [Ao,0] where Aq := —1; see Figure]^ in this sense, perturbations that reduce 
A are excitatory while those that increase A are inhibitory. The coupling and the choice of 
parameters will be made so as to construct a network in phase space where each connection 
goes once around a hysteresis loop within a subspace Pi (dehned in the following section). 


2.1 Dynamics of the model 

I I 

System (5) has symmetries Z2 ^ given by yk e-)■ —yk for each k and so the system is equivariant 
under the action of the group 




(k) 

2 • 


k=l 


We prove the existence of networks in phase space that realise the given graph and are robust 
to perturbations that respect this symmetry. To this end we denote by the subgroup of 
S corresponding to and dehne the following subspaces of phase space 


and 


:= £x(S£) = {(p, y) : = 0} 

■= nfc^£fix(Sfc) = {{p,y) : 2/fc = 0 

Pi ■= {{p, y) ■ 2/fc = 0 if fc 7^ £ and pj = 0 if j 7^ «(^) or u{i)} 
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Figure 6 : Bifurcation diagram of y = g{y, A) (given in ([^) for ?/ > 0 and A G M. 
there is a region of bistability between the pitchfork bifurcation at A = Aq := — 1, 
the saddle-node bifurcation at A = 0, ?/ = 1. 


Note that 
y = 0 and 


for £ = 1,..., Ue- The sets and Wg are invariant for all S-equivariant perturbations, while 
the Pi C Wi are invariant for the flow generated by system ([^ but not for arbitrary S- 
equivariant perturbations. However, we will show that there are connecting orbits in the Pi 
subspaces that are robust to small S-equivariant perturbations that preserve the invariance 
of the larger subspace Wi. We dehne 

5:={(p,2/) : = 1} ^ X (8) 


which is an n.y—1-dimensional sphere in the p-coordinates. This is invariant and normally 
attracting for F > 0, and so persists for appropriate choice of the parameters. We interpret 
the parameters A, B, C, D, F in (|^ as follows: 

• The constant A determines the default dynamics of the variables yk- we assume H > 0 
so that the equilibria are globally stable for F = C = 0. 

• The constant F > 0 determines how much is destabilised by there being a connection 
from that state. Let a{k) = j. Then if F > H -|- 1, will be linearly unstable in the 
yk direction. IfH-|-l>F>0 then the state will be linearly stable but excitable in 
the yk direction. 

• The constant C > 0 determines the mutual inhibition of the yk variables and suppresses 
more than one hysteresis loop becoming active at any time. 

• The constant F > 0 sets the rate of attraction to the equilibria in directions tangent 
to S. 


The constant F > 0 is set to lie within a range (relative to D) so that when one of the 
yk is active then there is a connection from ^a(k) to (see Figure]^ and equation (12) 
for details). 


• The constant F > 0 sets the rate of attraction of the p dynamics towards S. 
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Lemma 2.1 The system ^ has an invariant set S (defined in ([^j; for F > 0 this set 


attracts a neighbourhood of S. 

Proof: We show that if A := then A —)■ 1 as t —)■ cx) for typical initial 

conditions. Note that 






= FA(l-A). 

Hence on a timescale determined by F > 0, we typically have A —)■ 1 as t —)■ cxd. The only 


initial conditions where this is not the case will have pj = 0 for all j. 


QED 


2.2 Realisation of a graph as a heteroclinic network 

The following lemma shows that for an open region in parameter space the dynamics of 
system (|^ embeds the graph T as a heteroclinic network. 

Lemma 2.2 Consider the system 0 with equilibria at for j = 1,... ,n„. There is an 
open set of A, B, C, D, E, F such that for each £ = 1,..., rig there is a connecting orbit from 
ia(i) to within the three dimensional invariant subspace Pi. 

Proof: For ease of exposition, and without loss of generality, suppose i = I, a{i) = 1 and 
u{i) = 2. The system within the invariant subspace Pi can be written 


Pi = Pi{F{l - p^) + D{pIp^ - p'^)) - Eylpip2 
P 2 = P2{F{1 - p^) + D{pIp^ - p^)) + Eylp\ 

Vi = 9{yi,A-Bpl) 


(9) 


where p^ = p\+ pI.- Then the linearized stability within Pi is given by 

/ F{1 - 3pj - pI) + D{3pIp^ - pj) - Ep 2 yf pi{D{2plp2 - 4pj) - 2Fp2 - Eyf) 
pi{D{2pj - Ap 2 pI) - 2Fp2 + 2Eyl) F{1 - Spj - pj) + D{3plp^ - pf) 


yi(A - 2Bpi) 


0 



where g'(y, A) : 



dy 


0 0 
-D 0 
0 B-1-A 


-2F 

0 

0 
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while for ^2 = ( 0 , 1 , 0 ) it becomes 


-DO 0 
0 -2F 0 

0 0 - 1-/1 

Hence, we choose F>0, D>0, H>0 and D > 1 + H so that is a saddle with unstable 
direction (0, 0,1) and ^2 is a stable node. Observe that for this choice both and ^2 are 
hyperbolic and that all other eigenvalues in the direction of other pk are —D and therefore 
stable. 


As the subset C = SCiPi (where p\+p\ = 1) is attracting and invariant (by Lemma 2 . 1 ), 
we consider the dynamics on C parametrized by {9,yi) where pi = cos6, p 2 = sin^ and 
9 e [0,27r). From equation (|^ we have 


M 

dt 


^ = PiP 2 D{pI— pi) + Epi{pl+pl) = D sin 9 cos 9 (9) — cos^ (9)) + Eyf cos 9 

D , 

= — — sin46'+ D?/^cos 6 ' 


so that in this subspace we have 


d9 E> 2 n 

— = — sm 40 + Ey cos 9 
dt A ^ 

% = -yiiy"^ - if + A - B cos^ 9) 
dt 


( 10 ) 

( 11 ) 


where we drop the subscript from yi = y for notational convenience. 

We use (10) and ( [II| to deduce conditions on the parameters A, B, D, E that guarantee 
existence of a saddle-to-sink connection from ( 0 , y) = ( 0 , 0 ) to ( 0 , y) = (vr/ 2 , 0 ) corresponding 
to existence of a heteroclinic connection from ,^1 to .^2 within Pi. Note that the 0 = 0 nullclines 
are at cos 0 = 0 or at = D sin 0(2 cos^ 0 — 1) (shown by the dashed lines in Figure [^. For 
D > 0 the latter has a unique maximum in [ 0 , 7 r/ 2 ] at 


0 := sin"' ( — 


y ■= 


>D^/Q 

9E 


If D > 0 and D > 0, 0 > 0 whenever y > y. We require that y < 1 in order to rule out the 
possibility of any equilibria in y > 1 . 

The y = 0 nullclines are at y = 0 and B cos^ 0 = (1 — j/^)^ + A. If D > 1 + A the latter 
curve has a minimum in 0 at ?/ = 0 and maxima in 0 at ?/ = ±1 (dashed-dotted line in 
figure l^a)). Suppose that the line y = y hits the |/-nullcline first in [0, 7 r/ 2 ] at 0; this is given 
by 


cos 


0 = - 



+ A 
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Figure 7: The phase plane for {6,y) showing the nullclines in the region 6 G [0,7r/2]. 
(a) The dashed line shows 0 = 0 while the dash-dotted line shows yi = 0; equilibria are 
indicated with disks. If the nullclines have the given topology and are such that ^ < 1 and 
7 r/4 < 6 < 7r/2 then there will be a heteroclinic connection as shown (schematically as a 
solid red line). This can be achieved by choosing constants that satisfy (12); see text for 
more details, (b) For larger values of B the heteroclinic connection becomes an excitable 
connection (schematically as a solid red line) for amplitude 5 (schematically as a solid red 
line) and a new saddle equilibrium with y ^ 0 appears. 


If 7 r /4 < 0 < 7 r /2 and 0 < ^ < 1 then Figure [^a) shows that there will be a connection 
as desired; hence some sufficient (but by no means necessary) conditions for there to be a 
connection can be expressed as: 


0 < A <—, 1 + A<B, 0<Z1, 0<E and 



<D < 


9E 

Vo 


( 12 ) 


Note that C, E do not affect this argument; however C needs to be chosen to positive and 
large enough to avoid spurious connections to other stable dynamics and E needs to be 
positive for Lemma 2T to hold. Figure illustrates that this set is non-empty and open. 

QED 


To illustrate the effect of one component of the y dynamics becoming active and giving 
a connection, we show in Figure]^ the dynamics in pi and p 2 for (a) yi = 0 and (b) yi ~ 1. 
We summarise the construction above in the following Theorem: 


Theorem 2.3 Given any finite directed graph, there is a non-empty and open set of pa¬ 
rameter values A, B, C, D, E, E such that the system realises this graph as a heteroclinic 
network in a way that is robust to all perturbations to the eguations that preserved the sym¬ 
metries S. 


Proof: The previous calculations and Lemma 2.2 show that for the parameter region iden- 
tihed in ( 12 ) and with (F > 0 , F > 0 there are hyperbolic saddles that are connected 
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Figure 8: The parameter region that satishes the conditions (12) that permit a heteroclinic 
realization are illustrated here: (a) {A, B) can be chosen from the region shown in wavy 
lines, (b) For each choice there is an M < L such that D, E can be chosen from the region 
shown in wavy lines, where L = 9/a/6 and M = L(1 — ^jB/2 — A) which depends on A 
and B. The grey shaded region in (a) can be added to the allowable conditions if we permit 
excitable realizations with small 5 > 0. Note that not all of the networks are attracting, but 
for small enough B and suitable choices of C, F, numerical simulations indicate that they 
are. 



Figure 9: The dynamics of pi and p 2 on the unit circle in the case that yi is associated with 
a connection from di to t) 2 , i-e. in the case «(!) = 1, Ci;(l) = 2. (a) For y = 0 note that both 
and ^2 are stable nodes; note that is unstable in the yi direction, (b) For yi cs 1. As yi 
is increased from 0, the terms multiplied by E in ([^ removes (and all symmetric images) 
in a series of bifurcations. 
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by heteroclinic connections in the three-dimensional subspaces Pi for each These connec¬ 
tions are robust to perturbations that preserve the symmetry S because within the hxed 
point subspace Yi the equilibrium ^^(i) has one unstable direction and all other directions 
are stable while S,u}(e) is a sink. Moreover Pi <Z Yi and so there is a connection that is of 
saddle-sink type. Transversality of this connection means that it is robust to S-equivariant 
perturbations. QED 


We conjecture that, in the case of a strongly connected graph T, for some open subset 
of the set of parameters in Theorem 2.3 this heteroclinic network realises the graph as part 
of an asymptotically stable attractor that is a compact, chain recurrent invariant set. As in 
[ 0 ], the large “embedding attractor” will typically contain extra equilibria and connections 
but we conjecture that the proportion of time that typical trajectories visit equilibria that 
do not correspond to those in T will be very small and may go to zero as noise amplitude 
decreases; the larger attractor may be “invisible” [H] except on a subset that corresponds 
to an embedding of T. 

Some numerical evidence for these conjectures is given in the next section for a specific 
example. In brief justihcation, if has connections to .^2 and .^3 via yi and y 2 then for large 
enough C, almost every trajectory on is a connection to one of .^2 or ,^ 3 . This is 

suggested by the dynamics of j/i, y 2 which are governed by: 


^ = -y,({yl-lY + l + A-Bp\ + Cyl) (13) 

^ = -y,({yl-lf + l + A-Bp\ + Cy\). (14) 


Fixing pi = 1 and 1 + A — B = {] and examining the phase plane for this system, if C >2/3 
then all orbits are bounded in forwards time and the only attractors for this system are in 
1/1 = 0 and 1/2 = 0. This is preserved for 1 -|- A — i? close to zero. 


2.3 Realisation of a graph as an excitable network 

We give an additional result to show that a given graph can be realized as an excitable 
network for amplitude 5 > 0 for a range of parameter values: 

Theorem 2.4 Given any finite directed graph, there is an open set of parameter values 
A, B, C, D, E, F determining a minimum amplitude (threshold) 6th > 0 such that I® realises 
this graph as an excitable network with amplitude 6 for 6 > 5th but not for 5 < 5th- 

Proof: We choose A, B, C, D, E, F as in Theorem |2.3| except we do not require 1 + A — B < 0 
meaning that the yi nullcline may be detached from yi = 0; see Figure [^b). That is, the 
curve B cos^ 6 * = (1 — y‘^Y -|- A is undehned for 0 < // < h, where it is simple to show that 

5 = \l I- ^/B^; 
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Figure 10: (a) The Kirk-Silber network: five edges Cj connect four vertices t)*; (b) schematic 
of the coupling architecture of the nine-cell network realizing (a); the globally inhibitory 
connections are not shown. 

see the dot-dashed curve in Figure [^b). In this case is a sink within the invariant subspace 
Pi and there is a nearby saddle whose stable manifold forms part of the boundary of 
hF®(^i). Let ^0 be the smallest distance from to the stable manifold of Ci. 

More precisely, if we consider — = then we can estimate the 

closest approach of hF'^(Ci) to by 5, the point where the yi nullcline intersects 0 = 0. This 
gives 

^th = \j^ +0{y). (15) 

for small z/; in other words, for 6 > 6th there will be a connection for amplitude 6 while for 
6 < 5th we have i? 5 (^i) C and there is no connection for amplitude 5. QED 

In the case of multiple outgoing directions there will be multiple directions with a thresh¬ 
old of 5th and a similar argument to that following Theorem |2]^ suggests that for large enough 
C > 0 and 5 > 5th, almost all points in Bs{^k) are either in W^{^k) or in for some 

that is connected via one of these outgoing directions. 

3 Design of a system possessing a Kirk—Silber cycle 

One of the simplest examples of a network that shows competition between two heteroclinic 
cycles is the network of Kirk and Silber [18] where two order-three cycles (similar to that 
in Figure share a common edge. This is a useful system to understand how the system 
switches at the “decision point” in response e.g. to noise of differing amplitudes in different 
components. In a forthcoming paper [7] we explore the statistics of the switching process 
in terms of escape processes simultaneously along a number of heteroclinic or excitable 
connections; here we indicate some of the issues in this example. 

Let us consider a specific example of a realization of the network shown in Figure [T^a) 
using the coupled cell network illustrated in (b). To this end we consider the model perturbed 
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by additive noise wf 


pi = Pi{F{l 
p2 = P2{F{1 
p3 = P3{F{1 
Pi = Pa{F {1 
yi = 9{yi,A 

m = 5'(i/2, a 
m = 9{y3, A 
Vi = g{yi, A 
95 = 9 ( 95 , A 


-p ) + JJpPiP -p }} + J^[-y^piP2 + I/ 4 P 3 + 2 / 5 P 4 J - 

- + D{pIp^ - p^)) + E{-yjp2P3 - ylp 2 Pi + ylp\ 

- p^) + D{pIp‘^ - p^)) + E{-ylp3pi + p^p^) + ppW3 

- p2) + D{pIp‘^ - p^)) + E{-ylpiPi + ylpl) + p^Wi 

- Bpl + Ciy"^ - yl)) + r/iW5 

- Bpl + C{y‘^ - yl)) + r/ 2 W 6 

- Bpl + Cpy'^ -yl))+py,wj 

- Bpl + Ciy"^ - yl)) + r/4W8 

- Bpl + C{y^ - yl)) + r/5W9 


(16) 


where we choose the parameters as in (|^ except we allow different noise amplitudes in the p* 
directions. Figure 11 (a) and (b) show a time series for this case, where the parameters are 
chosen so that a heteroclinic network exists. Figure [^e ) shows a histogram of the residence 
times near equilibria for a much longer time series. 

Stone and colleagues [25l ES] have shown that for a heteroclinic cycle, the mean residence 
time near equilibria scales like (l/A)log(p) for small noise amplitude p, where A is the 
expanding eigenvalue at the equilibrium. We note that for the data shown in Figure [T^e), 
the residence times (not shown) near equilibrium ^2 are smaller than near the other three 
equilibria: this is due to the two possible ‘escape routes’ from that equilibrium. 


3.1 Bifurcation to an excitable Kirk-Silber network 

If B > A + 1, the equilibria E are connected to form an excitable network. We consider this 


same example (16) with parameters as in (j^ except for B = 1.49, and pj = 3 x 10 
this case there is an excitable connection with threshold 


-5 


In 


6th \Ji-VB-J ^ 0.07071 

from (15). Figure 11 (c) and (d) show a time series for this case. Figure [IT|(f) shows a 
histogram of the residence times near equilibria for a much longer time series. Similarly to 
in the heteroclinic case, we again note that the residence times near E appear to be smaller 
than near the other equilibria due to the presence of two escape routes. We also note that 
the shape of the distributions of escape times in the heteroclinic and excitable cases appears 
to be different. 

In particular, very long residence times are more likely in the excitable case than in 
the heteroclinic case (the distribution has a fatter tail). This can be seen in comparison of 
Figure [II|( e) and (f); although the timescales are different (due to the different parameters), 
we have scaled the a:-axes so that the mean residence times appear at the same position in 
each figure (these are approximately 10.3 and 42.3 respectively), and the fatter tail in (f) 
can be clearly seen, along with a much more peaked distribution in (e). 
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(b) 


(d) 



Figure 11: Time series and histograms of residence times for the Kirk-Silber example in (16). 
(a), (b) and (e) are for the heteroclinic case (parameters as in (|^, but with rjj = 3 x 10“® 
(j = 1,... ,5)). (c), (d) and (f) are for the excitable case (parameters as in (Q, but with 
B = 1.49 and rjj = 3x 10“^ {j = 1,, 5)). (a) and (c) show time series for the pj (linestyles 
are: pp. blue solid, p 2 . red dashed, ps magenta dotted, pp. black solid); (b) and (d) show 
time series for the yj, we only show j /2 (blue solid line) and (red dashed line), (e) and 
(f) show histograms of residence times near equilibria for each case for a much longer time 


series. 
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In contrast to the case for the heteroclinic network, we expect the residence times for 
the excitable network will be governed by a Kramers-type of escape process giving rise an 
exponential tail in the residence times. This and other statistical properties of switching 
and residence times near heteroclinic and excitable networks are being investigated for a 
forthcoming work |^. 


4 Discussion 


The main resnlts of the paper are the constructive model ([^ and Theorems 2.3 and 2.4 that 
give open sets of parameter values where one can robustly realize an arbitrary directed graph 
as an attracting heteroclinic network or an excitable network in phase space. 

Although the noise-free dynamics of the network is to some extent trivial (there will 
be slowing-down heteroclinic dynamics visiting a sequence of equilibria that depends on 
the initial condition for the heteroclinic network, or, the trajectory will remain at the hrst 
stable equilibrium for the excitable network) the dynamics of the networks become much 
more interesting in the presence of noise. It is known that addition of low noise to an 
attracting heteroclinic network can lead to random switching around a heteroclinic network 
in a temporally fairly regular manner [6] ; similarly, addition of noise to an excitable network 
leads also to random switching around an excitable network. We will explore elsewhere the 
switching probabilities and residence times near equilibria as a function of the noise strengths 
and parameter values [7j. 

Similarly, addition of very small inputs in the form of impulses to the yu variables allows 
one to control transitions between states in a way that depend on inputs and current state - 
and so perform hnite state computing in the system (|^ in a manner similar to [31 El 1231 EZj • It 
will be interesting to explore the computational potential of this network. Our construction 
is “wasteful” in the sense that only one cell will be active at any time; the encoding of states 
is very sparse compared to what nervous systems presumably achieve. It will be a challenge 
to see whether this construction can be adapted to achieve more dense encoding without 
losing the high level of control of the dynamics. 

Our concept of an excitable network (in phase space) needs to be distinguished from the 
more general concept of a network of excitable systems. A network of the latter type may or 
may not realise the former as a network in phase space, depending the nature and strength 
of the coupling. For example, networks of excitable systems (see for example [T3l fT6] ) may 
have many excitable states corresponding to various combinations of cells being active. 

The boundary between heteroclinic and excitable network dynamics for (|^ is on the line 
B = 1 + A and corresponds to a subcritical pitchfork bifurcation of the equilibria within 
the invariant subspace for each outgoing direction from By considering and 
(i.e. A, B depending on i) one can clearly design networks using (|^ that mix heteroclinic 
and excitable connections with thresholds 6i that may vary from one connection to another. 

The presence of microscopic noise in the heteroclinic network will result in trajectories 
wandering around the embedded graph with random choice of outgoing edges at each node. 
For A and B independent of k and low amplitude noise, this will appear to be a one-step 
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Markov process with a distribution of residence times. However, varying and in the 
heteroclinic network case introduces the possibility of “lift-off” and “memory” of the system 
trajectories in ([^ as discussed for a related system in [6]. If the transition past a node is 
sufficiently fast compared to previous nodes the transition probabilities between nodes may 
depend not just on the current state but on previously visited states; this will be dependent 
on the eigenvalues of the equilibria This gives the possibility of designing a system with 
more complex time correlations than a one-step Markov process. 
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